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Abstract 

We present new results building on the conservative deterministic spec- 
tral method for the space homogeneous Boltzmann equation developed by 
Gamba and Tharkabhushaman. This approach is a two-step process that 
acts on the weak form of the Boltzmann equation, and uses the machin- 
ery of the Fourier transform to reformulate the collisional integral into a 
weighted convolution in Fourier space. A constrained optimization prob- 
lem is solved to preserve the mass, momentum, and energy of the resulting 
distribution. 

Within this framework we have extended the formulation to the case of 
more general case of collision operators with anisotropic scattering mech- 
anisms, which requires a new formulation of the convolution weights. We 
also derive the grazing collisions limit for the method, and show that 
it is consistent with the Fokker-Planck-Landau equations as the grazing 
collisions parameter goes to zero. 

1 Introduction 

There are many difficulties associated with numerically solving the Boltzmann 
equation, most notably the dimensionality of the problem and the conservation 
of the collision invariants. For physically relevant three dimensional applica- 
tions the distribution function is seven dimensional and the velocity domain is 
unbounded. In addition, the collision operator is nonlinear and requires evalu- 
ation of a five dimensional integral at each point in phase space. The collision 
operator also locally conserves mass, momentum, and energy, and any approx- 
imation must maintains this property to ensure that macroscopic quantities 
evolve correctly. 
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Spectral methods are a deterministic approach that compute the collision 
operator to high accuracy by exploiting its Fourier structure. These methods 
grew from the analytical works of Bobylev [2] developed for the Boltzmann 
equation with Maxwell type potential interactions and integrable angular cross 
section, where the corresponding Fourier transformed equation has a closed 
form. Spectral approximations for these type of models where first proposed 
by Pareschi and Perthamc [11]. Later Pareschi and Russo [12] applied this 
work to variable hard potentials by periodizing the problem and its solution 
and implementing spectral collocation methods. 

These methods require 0(N 2d ) operations per evaluation of the collision 
operator, where iV is the total number of velocity grid points in each dimension. 
While convolutions can generally be computed in 0(N d log N) operations, the 
presence of the convolution weights requires the full 0(N 2d ) computation of 
the convolution, except for a few special cases, e.g., the Fokker-Planck-Landau 
collision operator [13, 10]. Spectral methods provide many advantages over 
Direct Simulation Monte Carlo Methods (DSMC) because they are more suited 
to time dependent problems, low Mach number flows, high mean velocity flows, 
and flows that are away from equilibrium. In addition, deterministic methods 
avoid the statistical fluctuations that are typical of particle based methods. 

Inspired by the work of Ibragimov and Rjasanow [8], Gamba and Tharkab- 
hushanam [6, 7] observed that the Fourier transformed collision operator takes 
a simple form of a weighted convolution and developed a spectral method based 
on the weak form of the Boltzmann equation that provides a general framework 
for computing both elastic and inelastic collisions. Macroscopic conservation 
is enforced by solving a numerical constrained optimization problem that finds 
the closest distribution function in L2 to the output of the collision term that 
conserves the macroscopic quantities. 

This paper presents an extension of this method to elastic collisional models 
that arc anisotropic with respect to the scattering direction. When the angular 
cross section of the collision kernel has a non-integrable singularity at zero, un- 
der certain conditions one can formally calculate the Landau asymptotics [9, 14] 
to obtain the well known Fokker-Planck-Landau equation. In particular, this 
can apply to the well known Rutherford Coulombic potential [15], which is of 
interest for plasma physics applications. These facts combined with our spectral 
constrained optimization computational schemes allow us to investigate a gen- 
eral class of soft potential collisional kernels with anisotropic singular angular 
sections as well as the grazing collision transition regime between the classical 
Boltzmann equation to the Fokker-Planck-Landau equation. 

2 The space homogeneous Boltzmann equation 

The space homogeneous Boltzmann equation is given by the initial value prob- 
lem 

|/(v,*) = ^W,/), (1) 



2 



with 



ver, /(«,0) = / (v) 



where /(v, t) is a probability density distribution in v-space and /o is assumed 
to be locally integrable with respect to v. The dimensionless parameter Kn > 
is the scaled Knudsen number, which is defined as the ratio between the mean 
free path between collisions and a reference macroscopic length scale. 

The collision operator Q(f,f) is a bilinear integral form in (v, t) given by 

Q(/,/)(v,t)= / / B(\v-v*\,co S 6)(f(v:)f(v')-f(v*)f(v))d<Tdv*, 

(2) 

where the velocities v',v' + are determined through a given collision rule (3), 
depending on v, v*. The positive term of the integral in (2) evaluates / in the 
pre-collisional velocities that can result in a post-collisional velocity the direc- 
tion v. The collision kernel i?(|v — v*|,cos#) is a given non-negative function 
depending on the size of the relative velocity u := v — v* and cos 8 — 3 r£, 
where a in the n — 1 dimensional sphere is referred ro as the scattering 

direction which also coincides with the direction of the post-collisional elastic 
relative velocity. 

For the following we will use the elastic (or reversible) interaction law 

v'=V+i(M<7-u), <=V*-i(|u|<7-u) (3) 

B(\u\, cosd) = \u\ x b(cos6) . 

The angular cross section function b(cos9) may or may no be integrable with 
respect to a on S d ~ x . When integrability holds is referred to as the Grad cut-off 
assumption on the angular cross section. 

The parameter A regulates the collision frequency as a function of the relative 
speed |u|. This corresponds to the interparticle potentials used in the derivation 
of the collisional kernel and are referred to as variable hard potentials (VHP) for 
< A < 1, hard spheres (HS) for A = 1, Maxwell molecules (MM) for A = 0, and 
variable soft potentials (VSP) for — 3 < A < 0. The A = — 3 case corresponds to 
a Coulombic interaction potential between particles. If 6(cos 9) is independent 
of g we call the interactions isotropic, e.g, in the case of hard spheres in three 
dimensions. 



2.1 Spectral formulation 

The key step our formulation of the spectral numerical method is the use of 
the weak form of the Boltzmann collision operator. For a suitably smooth test 
function </>(v) the weak form of the collision integral is given by 

/ 0(/,/)^(v)dv= / /(v)/(v.)B(|u|,cose)(^(v , )-^(v))dtrdv.dv 

JR d JR d xR d xS d - 1 

(4) 
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If one chooses 

0(v) = e-^/i^f, 
then (4) is the Fourier transform of the collision integral with Fourier variable 





Q{Q = —±— / Q(f,f)e-«-*dv 



/* « w, ..B(|u|,cos0) _ iC . v - -if-vNj j . 

= / /(v)/(v«) (e 4 -e * )d<rdv*dv 

= / G(u,CW(v)/(v-u)](C)du, (5) 

jR d 

where [•] = .F(-) denotes the Fourier transform and 

G(u,() = H x J 6(cos0)(e-^ c > |CT e^ c - u -l)rfa (6) 

Further simplification can be made by writing the Fourier transform inside the 
integral as a convolution of Fourier transforms: 

Q(C)= / G(£,C)/(C-0/m. (7) 

jR d 

where the convolution weights G(£, Q are given by 

%,C) = — =L- / G(u,C)e-^rfu (8) 
= _L / | M | A e- l « u / 6(cos6») f e -^Ok e ^C-u _ A d(jdu 

These convolution weights can be precomputcd once to high accuracy and stored 
for future use. For many collisional models, such as isotropic collisions, the 
complexity of the integrals in the weight functions can be reduced dramatically 
through analytical techniques. However in this paper we make no assumption 
on the isotropy of b and derive a more general formula. 

We begin with G(u, () and define a spherical coordinate system for a with 
a pole in the direction of u and obtain 

G(u,C) = 2tt\u\ x J\(cos6) S m0 (e'la—W-Jo ^M 8 ^^ _ ^ d0j 

(9) 

where Jo is the Bessel function of the first kind, j, k are unit length basis vectors 
that are mutually orthogonal to u, and C, 1 - = ( — (( ■ u/|u|)u/|u|. Note that for 
the isotropic case b(cos0) drops out and one can take £ as the polar direction 
for a, resulting in an explicit expression involving a sine function. 
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G is the Fourier transform of G, however we must take this transform on a 
ball centered at rather than simply taking the Fast Fourier Transform (FFT) 
[5] of G, which ensures that the weights are real- valued. 

The convolution weights G(C,£) in 3 dimensions are computed as follows 



G(C,£)=2tt 



= 2tt 



L 



b(cos9) sin# 
b(cos9) sin 9 



; if.u(i-co S <>) Jo fl\ u \\^\ sin e\ _ i 



dOdu 



e - l r( 4 -i(l-cose)).r, Jo Q r | C ±| sin ^ _ e - lr t- V 



dOd-qdr. 



We now take £ to be the polar direction for the spherical integration of r\ 
and use that G is real- valued to obtain 



G(C,£)=47r 2 



_A+2 



Jo 



b(cos9) sin sin </>J ( 



o r 



5 icr 



sin</> x 



cos I r(£ — ^(1 — cos 6*)) • |^|- cos^ Jo ^ r KI sin0sin# 



cos [ r£ ■ j^j- cos i 
(10) 



d9d(j)dr. 



This requires a three dimensional integral for the N 6 pairs (£,£)> which is 
two orders of magnitude more than the isotropic case, but as before this weight 
is precomputed only once and used in any subsequent computations with the 
collisional model. 



3 The Conservative Numerical Method 

3.1 Velocity space discretization 

In order to compute the Boltzmann equation we must work on a bounded veloc- 
ity space, rather than all of M. d . However typical distributions are supported on 
the entire domain, for example the Maxwellian equilibrium distribution. Even 
if one begins with a compactly supported initial distribution, each evaluation of 
the collision operator spreads the support by a factor of y/2, thus we must use a 
working definition of an effective support of size R for the distribution function. 
Bobylev and Rjasanow [3] suggested using the temperature of the distribution 
function, which typically decreases as exp(— \v\ 2 /2T) for large \v\, and used a 
rough estimate of R m 2y/2T to determine the cutoff. Thus, we assume that 
the distribution function is negligible outside of a ball 

B Rx (V(x)) = {v e R d : |v - V(x)| < R x }, (11) 

where V(x) is the local flow velocity which depends in the spatial variable x. 
For ease of notation in the following we will work with a ball centered at and 
choose a length R large enough that Br x (V(x)) C -Br(O) for all x. 
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With this assumed support for the distribution /, the integrals in (7) will 
only be nonzero for u £ B2r(0). Therefore, we set L = 2R and define the cube 

C L = {veR d : \ VJ \ <L, j = l,...,d} (12) 

to be the domain of computation. With this domain the comptuation of the 
weight function integral (10) is cut off at tq = L. 

Let N € N be the number of points in velocity space in each dimension. Then 
we establish a uniform velocity mesh with Aw = ^ and due to the formulation 
of the discrete Fourier transform the corresponding uniform Fourier space mesh 
size is given by A£ = J. 

3.2 Collision step discretization 

To simplify notation we will use one index to denote multidimensional sums 
with respect to an index vector m 

N-l N-l 

E = E • 

m— mi, ...,m(2=0 

To compute Q(Ck), we first compute the Fourier transform integral giving 
/(Cfe) via the FFT. The weighted convolution integral is approximated using the 
trapezoidal rule 

N-l 

Q(Ck) = E Ck)/(U)/(Ck - Cm)^ m , (13) 

m=0 

where w m is the quadrature weight and we set /(£k — £ m ) = if (£k — £ m ) is 
outside of the domain of integration. We then use the inverse FFT on Q to 
calculate the integral returning the result to velocity space. 

Note that in this formulation the distribution function is not periodized, as 
is done in the collocation approach of Pareschi and Russo [12]. This is reflected 
in the omission of Fourier terms outside of the Fourier domain. All integrals 
are computed directly only using the FFT as a tool for faster computation and 
the convolution integral is accurate to at least the order of the quadrature. The 
calculations below use the trapezoid rule, but in principle Simpson's rule or 
some other uniform grid quadrature can be used. However, it is known that the 
trapezoid rule is spectrally accurate for periodic functions on periodic domains 
(which is the basis of spectral accuracy for the FFT), and the same arguments 
can apply to functions with sufficient decay at the integration boundaries [1]. 
These accuracy considerations will be investigated in future work. The overall 
cost of this step is 0{N 2d ). 

3.3 Discrete conservation enforcement 

This implementation of the collision mechanism does not conserve all of the 
quantities of the collision operator. To correct this, we formulate these conser- 
vation properties as a Lagrange multiplier problem. Depending on the type of 



G 



collisions we can change this constraint set (for example, inelastic collisions do 
not preserve energy), but we will focus on the case of elastic collisions, which 
preserve mass, momentum, and energy. 

Let M — N d be the total number of grid points, let Q = (Qi, • • • , Qm) T be 
the result of the spectral formulation from the previous section, written in vector 
form, and let ojj be the quadrature weights over the domain in this ordering. 
Define the integration matrix 




where v l , i = 1, 2, 3 refers to the ith component of the velocity vector. Using this 
notation, the conservation method can be written as a constrained optimization 
problem. 



Find Q = (Qi, . . . , Q M ) T that minimizes -|| Q - such that CQ = (14) 
Formulating this as a Lagrange multiplier problem, we define 

M 
3=1 

The solution is given by 

Q = Q + C(CC T ) _1 CQ 

= PnQ (16) 
Overall the collision step in semi-discrete form is given by 

g = P,Q (17) 

The overall cost of the conservation portion of the algorithm is a 0(N d ) 
matrix- vector multiply, significantly less than the computation of the weighted 
convolution. 



3.4 Grazing collisions limit 

The Fokker-Planck-Landau equation is used to describe binary collisions occur- 
ring in a plasma, for example, and can be shown to be an approximation of 
the Boltzmann equation in the case where the dominant collision mechanism is 
that of grazing collisions, i.e., collisions that only result in very small deflections 
of particle trajectories, as is the case for Coulomb potentials with Rutherford 
scattering [15] between charged particles. The equation is given by 

Q L (f,f)=V v -( [ \ u f+\l-?^)(f(v*)V v f(v)-f(v)(V v f)(v*))dv?l . 

(18) 
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Applying the same methodology used to derive the spectral method for the 
Boltzmann equation, the Fourier transform of the Fokker-Planck-Landau equa- 
tion takes the form (5), where the weight function G(u, £) in (6) is now given 
by 

g(u,c) = 4H^Wu.c)-Jh 2 |C ± I 2 ). (19) 

To show that the spectral method for Boltzmann operator is consistent with 
this form of Fokker-Planck-Landau operator, we must take the grazing collisions 
limit, which requires that the angular scattering function is consistent with the 
singular rates of Rutherford scattering. Indeed, it is enough to assume that the 
collision kernel satisfies the following. 

Let e > be the small parameter associated with the grazing collisions limit. 
A family of kernels b e {9) represents grazing collisions if: 

• lim e ^ A e = 2-7T b(cos9, s)(l — cos 9) sm9d9 = A < oo 

• lim e ^ = 2tt r o 7r 6(cos6»,e)(l -cos9) 2+k sm 9 d9 -> for fc>0. 

In particular, the angular part of the classical Rutherford scattering for 9 
near corresponds to a family b s (cos9) given by 

&£ 

b(cos0,e)sm9 = — -1 0>£ , (20) 

which satisfies the above conditions for A = 8. This is the potential we take for 
the following calculation and numerical simulation of the grazing collision limit 
for the Boltzmann equation for Coulombic interactions in three dimensions, i.e., 
A = -3. 

For the grazing collisions limit, take Taylor expansion of the exponential 
term in G (9) 

e i|((l-co S e)C.U+|«|C.jsine S in0+|„|C-fesin0co S 0) _ j = _ cog . u + ^ . j s [ n Q sill(j) + \u\( ■ k sin 9 COS ( 

«2 

- ^-((1 - cos 0)C • u + |u|C • j sin 9 sin <j> + |«|C ■ k sin 9 cos (f>) 2 + (9(|C| 3 \u\ 3 9 3 ) 
8 

We can toss the 9 3 terms because as e — > 

6 e (cos6») sin 9 9 3 d9^ 

'o 

Taking the </> integral of the first order term gives 

in (3(1 — cos9)( ■ u. 

For the second order term, taking the <j> integral results in, after removing inte- 
grals that evaluate to zero by symmetry, 

/'"(I 1 - C0S ^) 2 (C ' u) 2 + |u| 2 (C • jf sin 2 9 sin 2 <\> + \u\ 2 (( ■ kf sin 2 9 cos 2 </>) 



P 71 

Jo 
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We can neglect the (1 — cos#) 2 term as it is 0(9 4 ). Thus we obtain the small e 
approximation for G 



H 2 IC 



2|^_L|2^2 

- sin 2 i 



G(u,C) = n \ u \ X J q b e (cos9)sm6 ^/3(1 - cos<9)C u 

From the earlier definition of the collision kernel, we have that 
27rsin#(l — cos9)b E (cos9) — > A <5e = o 

Thus 

G(u , = ^(tf(C.u)-H^!) (21 , 

When A = 8, this gives the same G(u, Q as found above by applying the 
spectral technique to the Fokker-Planck-Landau equation. 

3.5 Computing G for singular scattering kernels 

Numerically calculating the weights G to high accuracy can be difficult for 
singular scattering kernels, due to the precise nature of the cancellation at the 
left endpoint of the integral. Using b E (cos9) from (20), the 9 integral in (10) is 
given by 

-^-(cos(ci(l - cos6>) - c 3 )J (c 2 sin #) - cos(c 3 ))d9, (22) 

where C\,ci,c 3 depend on the current values of r, £, £ following from the full 
formulation of G. When e << 1 the bulk of the integration mass occurs near the 
left endpoint of the 9 interval, however this presents a challenge for a numerical 
quadrature package to compute. For 9 « 1 there is a subtraction of two nearly 
equal numbers, which causes floating point errors, as well as a multiplication 
of two numbers of very different magnitudes (e and 9~ A ). To alleviate this, we 
split the integration interval into two pieces, and use the first term of the Taylor 
expansion of the integrand for 9 « 1: 

J (~~[' C0S ( C 3) + y sin ( c 3)^ d9+ J -^(cos(ci(l-cos^)-c 3 )Jo(c2sin^)-cos(c 3 ))^. 

(23) 

These integrals are computed using the GNU Scientific Laboratory integration 
routines [4]. We use cquad to compute the first 9 integral, which has proved 
to be most stable choice for this near-singular integrand. The adaptive Gauss- 
Konrod quadrature qag is used for all other integrals used in computing the 
weights G. 

4 Numerical results 

To illustrate that this method captures the correct behavior for grazing colli- 
sions, we take a Coulombic potential (A = —3) and set e = 10~ 4 . Similar to 
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what was done in [13] for the Landau equation, we take the axially symmetric 
initial condition 



1.5 



/(v,0) = 0.01exp^-10 ^ l5 

We take a domain size of L = 5 and compute to time t = 100 with a timestep 
of 0.001. The results are shown in Figure 4. 




Figure 1: Slice of distribution function at different times. The solid dots are 
the values of the Maxwellian distribution associated with the initial data. The 
solid lines are a spline reconstruction based on the grid values (hollow circles) . 



5 Conclusions and future work 

We have derived the precomputed weight formulas for the more generalized 
anisotropic collisional models within the Boltzmann equation. We also showed 
that the spectral method for the Boltzmann equation is consistent with the lim- 
iting Fokker-Planck-Landau equation under suitable assumptions on the scat- 
tering kernel. One other important note is that this method could be a good 
candidate for collisional models where the collision mechanism is unknown and 
only experimentally determined. In addition, as the Landau equation is used to 
model collisions of charged particles in plasma we will seek to add field effects 
to the space inhomogeneous Boltzmann equation, resulting in the Boltzmann- 
Poisson or Boltzmann-Maxwell systems. The inhomogeneous method uses op- 
erator splitting between the collision and the transport terms, so in principle 
one can use whatever transport solver one wants for the spatial terms in the 
equation. 
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